################################################################
### Replication file for the MC experiments for "Long-Term Effects for Models with Temporal Dependence"
###
### Code created by dave carter & curt signorino, 4/27/2007 and modified by LKW
################################################################


source("Monte Carlo Programs.R")


################################################################
### Common parameters
################################################################
maxperiods=140
t=1:(maxperiods+10)
       
b1 = 1

xrange = c(-2, 2)

sims <- 1000
no.models <- 6
maxt <- 40
kn <- c(1, 4, 7)


################################################################
### Negative hazard rate (based on a Weibull, B-S Jones p. 25)
################################################################
set.seed(649)

lam = 1
rho = 0.25
h = lam*(rho)*(lam*t)^(rho-1)
h = h[2:length(h)]*12

### MC Scenarios: N and b.0 (to control for % of 1s in Y)
S <- matrix(c(1000, 1000, 5000, 5000, 10000, 10000,
		  -3, -2, -3, -2, -3, -2), byrow = TRUE, nrow = 2)

for(s in 1:ncol(S)){
   Nobs = S[1,s]
   b0 = S[2,s]

   print(paste0("Negative: Scenario #", s, " of ", ncol(S)))

   ### Simulate the data
   simtd("Negative/Negative", sims, no.models, maxt, kn, sc = s)

   ### Analyze the Monte Carlo experiments
   analtd("Negative/Negative", b0, b1, sc = s)

   ### Analyze the hazard rate
   analh("Negative/Negative", Nobs, b0, b1, maxt, sc = s)

   ### Examine first differences of 1-unit shift in X: average predictive differences
   analapd("Negative/Negative", Nobs, b0, b1, maxt, sc = s)

}

################################################################
### Positive hazard rate (based on a Weibull, B-S Jones p. 25)
################################################################
set.seed(648)

lam = 0.15
rho = 2.5
h = lam*(rho)*(lam*t)^(rho-1)

### MC Scenarios: N and b.0 (to control for % of 1s in Y)
S <- matrix(c(1000, 1000, 5000, 5000, 10000, 10000,
		  -3, -2, -3, -2, -3, -2), byrow = TRUE, nrow = 2)

for(s in 1:ncol(S)){
   Nobs = S[1,s]
   b0 = S[2,s]

   print(paste0("Positive: Scenario #", s, " of ", ncol(S)))

   ### Simulate the data
   simtd("Positive/Positive", sims, no.models, maxt, kn, sc = s)

   ### Analyze the Monte Carlo experiments
   analtd("Positive/Positive", b0, b1, sc = s)

   ### Analyze the hazard rate
   analh("Positive/Positive", Nobs, b0, b1, maxt, sc = s)

   ### Examine first differences of 1-unit shift in X: average predictive differences
   analapd("Positive/Positive", Nobs, b0, b1, maxt, sc = s)
}




################################################################
### Non-monotonic parabolic failure rate (formula from CS R script file)
################################################################
set.seed(648)

lam = 0.25
h = 1/3*(lam/2*(t-16))^2
h = h*2

### MC Scenarios: N and b.0 (to control for % of 1s in Y)
S <- matrix(c(1000, 1000, 5000, 5000, 10000, 10000,
		  -3, -2, -3, -2, -3, -2), byrow = TRUE, nrow = 2)

for(s in 1:ncol(S)){
   Nobs = S[1,s]
   b0 = S[2,s]

   print(paste0("Non-Monotonic 1: Scenario #", s, " of ", ncol(S)))

   ### Simulate the data
   simtd("Non-Monotonic 1/NM1", sims, no.models, maxt, kn, sc = s)

   ### Analyze the Monte Carlo experiments
   analtd("Non-Monotonic 1/NM1", b0, b1, sc = s)

   ### Analyze the hazard rate
   analh("Non-Monotonic 1/NM1", Nobs, b0, b1, maxt, sc = s)

   ### Examine first differences of 1-unit shift in X: average predictive differences
   analapd("Non-Monotonic 1/NM1", Nobs, b0, b1, maxt, sc = s)
}


################################################################
### Non-monotonic (log-logistic)
################################################################
set.seed(648)

lam = 0.25
rho = 1.5
h = (lam*rho*(lam*t)^(rho-1)) / (1+(lam*t)^rho)
h = h*10

### MC Scenarios: N and b.0 (to control for % of 1s in Y)
S <- matrix(c(1000, 1000, 5000, 5000, 10000, 10000,
		  -3, -2, -3, -2, -3, -2), byrow = TRUE, nrow = 2)

for(s in 1:ncol(S)){
   Nobs = S[1,s]
   b0 = S[2,s]

   print(paste0("Non-Monotonic 2: Scenario #", s, " of ", ncol(S)))

   ### Simulate the data
   simtd("Non-Monotonic 2/NM2", sims, no.models, maxt, kn, sc = s)

   ### Analyze the Monte Carlo experiments
   analtd("Non-Monotonic 2/NM2", b0, b1, sc = s)

   ### Analyze the hazard rate
   analh("Non-Monotonic 2/NM2", Nobs, b0, b1, maxt, sc = s)

   ### Examine first differences of 1-unit shift in X: average predictive differences
   analapd("Non-Monotonic 2/NM2", Nobs, b0, b1, maxt, sc = s)
}



################################################################
### Negative hazard rate (based on a Weibull, B-S Jones p. 25)
### Non-proportional hazard: positive (tt/nph)
################################################################
set.seed(648)

maxperiods=140
t=1:(maxperiods+10)
       
b1 = 1

detach(tddata)

xrange = c(-2, 2)

sims <- 1000
no.models <- 6
maxt <- 40
kn <- c(1, 4, 7)

Nobs = 1000
b0 = -3

lam = 1
rho = 0.25
h = lam*(rho)*(lam*t)^(rho-1)
h = h[2:length(h)]*12
s = 1

### MC Scenarios: N and b.0 (to control for % of 1s in Y)
S <- matrix(c(5, 10, 25, 50,
		  -3, -3, -3, -3), byrow = TRUE, nrow = 2)

detach(tddata)

for(s in 1:ncol(S)){
   nphv = S[1,s]
   b0 = S[2,s]

   print(paste0("Negative with Positive NPH: Scenario #", s, " of ", ncol(S)))

   simtd.nph("Negative/NPH Positive/NPH Positive", sims, no.models, maxt, kn, sc = s, nph = nphv)

   ### Analyze the Monte Carlo experiments
   analtd("Negative/NPH Positive/NPH Positive", b0, b1, sc = s)

   ### Analyze the hazard rate
   analh("Negative/NPH Positive/NPH Positive", Nobs, b0, b1, maxt, sc = s)

   ### Examine first differences of 1-unit shift in X: average predictive differences
   analapd.nph("Negative/NPH Positive/NPH Positive", Nobs, b0, b1, maxt, sc = s)
}


################################################################
### Positive hazard rate (based on a Weibull, B-S Jones p. 25)
### Non-proportional hazard: negative (tt/nph)
################################################################
set.seed(648)

maxperiods=140
t=1:(maxperiods+10)
       
b1 = 1

xrange = c(-2, 2)

sims <- 1000
no.models <- 6
maxt <- 40
kn <- c(1, 4, 7)

Nobs = 1000
b0 = -3

lam = 0.15
rho = 2.5
h = lam*(rho)*(lam*t)^(rho-1)

### MC Scenarios: N and b.0 (to control for % of 1s in Y)
S <- matrix(c(-5, -10, -25, -50,
		  -3, -3, -3, -3), byrow = TRUE, nrow = 2)

detach(tddata)

for(s in 1:ncol(S)){
   nphv = S[1,s]
   b0 = S[2,s]

   print(paste0("Positive with Negative NPH: Scenario #", s, " of ", ncol(S)))

   simtd.nph("Positive/NPH Negative/NPH Negative", sims, no.models, maxt, kn, sc = s, nph = nphv)

   ### Analyze the Monte Carlo experiments
   analtd("Positive/NPH Negative/NPH Negative", b0, b1, sc = s)

   ### Analyze the hazard rate
   analh("Positive/NPH Negative/NPH Negative", Nobs, b0, b1, maxt, sc = s)

   ### Examine first differences of 1-unit shift in X: average predictive differences
   analapd.nph("Positive/NPH Negative/NPH Negative", Nobs, b0, b1, maxt, sc = s)
}



